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ABSTRACT 

Radio telescopes are used to accurately measure the time of arrival (ToA) of radio 
pulses in pulsar timing experiments that target mostly millisecond pulsars (MSPs) due 
to their high rotational stability. This allows for detailed study of MSPs and forms the 
basis of experiments to detect gravitational waves. Apart from intrinsic and propaga¬ 
tion effects, such as pulse-to-pulse jitter and dispersion variations in the interstellar 
medium, timing precision is limited in part by the following: polarization purity of the 
telescope’s orthogonally polarized receptors, the signal-to-noise ratio (S/N) of the pul¬ 
sar profile, and the polarization fidelity of the system. Using simulations, we present 
how fundamental limitations in recovering the true polarization reduce the precision of 
ToA measurements. Any real system will respond differently to each source observed 
depending on the unique pulsar polarization profile. Using the profiles of known MSPs 
we quantify the limits of observing system specifications that yield satisfactory ToA 
measurements, and we place a practical design limit beyond which improvement of 
the system results in diminishing returns. Our aim is to justify limits for the front-end 
polarization characteristics of next generation radio telescopes, leading to the Square 
Kilometre Array (SKA). 

Key words: instrumentation: polarimeters — (stars:) pulsars: general — radio con¬ 
tinuum: general — techniques: polarimeters 


1 INTRODUCTION 

Any dual-polarization polarimeter is characterized by a de¬ 
gree of polarization purity, i.e. the cross-polarization be¬ 
tween orthogonal feeds, and the extent to which calibration 
can be used to retrieve accurate polarization information. 
In this paper we examine how both of these limitations af¬ 
fect pulsar time of arrival (ToA) measurements, especially 
in the case of millisecond pulsars (MSPs). We do this by 
simulating ToA measurements through a sampling of the 
signal-to-noise ratio (S/N), calibration error, and intrinsic 
polarization leakage parameter space. Intrinsic polarization 
leakage includes the apparent leakage between orthogonal 
receptors due to differential receptor gains, in addition to the 
cross-coupling between receptors, which is typically thought 
of as ‘pola-rization leakage’. 

Simulations are used, as it is difficult to analytically 
quantify the effects of calibration error, integration time, 


and intrinsic polarization leakage in a general form. The 
ToA measurement error depends on the ability to observe 
pulsar profiles with high fidelity. By profile, we mean the 
stable average shape of the radio pulse of a given pulsar, and 
its polarization properties. We perform our analysis using 
profiles from the 20 MSPs in Manchester et al. (20131. 

A fundamental limit to any ToA measurement is the de¬ 
sign of the polarimeter feeds, which is set by the telescope 
specifications. As pulsar timing is a key science project for 
Square Kilometre Array (SKA), see Janssen et al. (20151 
and Cordes et al. (20041, it is important to consider the sci¬ 
ence limitations set during the design process. The decime¬ 
tre wavelength band, where pulsars are typically observed 
for timing, will be covered by both dishes and aperture ar¬ 
rays. The analysis presented here applies to both telescope 
types. 

For the design of the feeds, we need to consider the ca¬ 
pacity of any dish or aperture array to produce data from 
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which the true polarization of a signal can be recovered. For 
example, the full polarization description of an incoming 
signal can never be recovered by a single dipole, no mat¬ 
ter how good the calibration procedure is. On more realistic 
systems, there is a fundamental limitation in recovering the 
polarization state due to differential gains between orthogo¬ 
nal receptors. These differential gains coupled with noise in 
the system result in measurement errors which can not be 
corrected via any currently used calibration procedure. This 
affects all To A measurement methods. For total intensity 
(e.g. Taylor (1992l) timing methods, i.e. the determination 
of a ToA through cross-correlation of a total power template, 
this is an effect in addition to calibration error. Techniques 


such as the invariant interval (Britton 20001 and matrix tem¬ 


plate matching (van Straten 


20061 methods, despite being 


largely independent of polarization calibration error, are also 
affected. 

The remainder of this introduction covers the relevant 
Jones and Mueller mathematical formalisms necessary to de¬ 
scribe intrinsic polarization leakage. In Sectionj^we describe 
the simulation setup and the strategy for exploring the rel¬ 
evant parameter space, how the simulated observations are 
generated, and the methods for determining the ToA. Re¬ 
sults are presented in Section Discussion of the results 
and the implications for current and future telescopes are 
presented in Section 


1.1 Intrinsic polarization leakage in Jones and 
Mueller formalism 

We are interested in describing the intrinsic polarization 
leakage of a linear feed, dual-polarization system as this is a 
typical design for single pixel feeds and phased arrays such 
as phased-array feeds (PAFs) and aperture arrays (AAs) 
used for pulsar timing. Intrinsic polarization leakage can 
be described with the mathematical structure developed for 
the radio interferometer measurement equation (RIME) pre¬ 
sented in Hamaker et al. (19961 and Smirnov ( 2011a| . Addi¬ 
tionally, the RIME can be extended to phased arrays where, 
to hrst order, the formed beam is a linear combination of the 
individnal element beams; a full description would also in¬ 
clude element mutual coupling terms. Jones matrix formal¬ 
ism is useful to frame the RIME in terms of instrumenta¬ 
tion and environmental effects on an electromagnetic signal. 
The Mueller matrix formalism, which is used in our simula¬ 
tions, is useful in interpreting the RIME in terms of detected 
power, represented by the Stokes parameters of the signal. 

In Jones formalism, transformations are applied to an 
input electromagnetic signal to produce the observed signal. 
The transformation from the complex electromagnetic sky 
Jones vector e to the observed voltage Jones vector v is 


V = Ja 


( 1 ) 


where Jays is the total system Jones matrix representation 
which is constructed out of multiple linear Jones transfor¬ 
mations, each of which can have dependence on time, ob¬ 
serving frequency, and source direction ( Smirnov|2011b |. In 
the scope of this paper we are interested in the effect of in¬ 
trinsic polarization leakage. The polarization leakage matrix 
D is usually defined as a direction-independent Jones ma¬ 
trix with the direction-dependent polarization leakage com¬ 
ponents incorporated into the primary beam matrix E. For 


this work we are not focusing on the primary beam direction- 
dependent sensitivity, but only the potentially direction- 
dependent polarization leakage, thus for phased arrays we 
are defining D to also include the direction-dependent po¬ 
larization leakage 

where (6, (f>) are reference frame dependent position angles. 
In Equation dp^q is the intrinsic leakage of feed p into 
feed q. In feed design, the off-diagonal terms are minimized. 
Ideally they are zero. For phased arrays, however, projection 
effects will cause D to vary with observing direction. 

We can then define an explicit RIME for our dish and 
phased array systems as 


d. 


'P^q (^5 ^5 0 ) 
1 


( 2 ) 


Jsys.dish = B(!/)G(f)CD( 


(3a) 


n 

Jsys,PA = Ew i{y, 9, 9, (j>) (3b) 

i = l 

where B is the frequency-dependent bandpass structure, and 
G is the time-dependent electronic gain. B and G are diago¬ 
nal matrices. The idealized nominal feed configuration C is 
a coordinate transform from the sky to observing frame. For 
a single pixel feed, on axis observation, D has no direction 
dependence and can be simplified to D(iz). 

For a phased array, the Jays is the weighted sum over all 
n elements of the array, where Wi are the complex weights 
to shape the array beam pattern to ‘point’ in a direction. 
For a PAF, the direction dependence of D will be relative to 
the dish pointing centre. For an aperture array, the direction 
dependence of D will be relative to the direction of boresight, 
usually zenith. 

Jones formalism is useful for understanding the instru¬ 
mental effects on a signal. For our simulation, however, we 
use pulsar profiles described as Stokes vectors, for which 
Mueller matrices are used to perform operations. Using the 
notation from Hamaker et al. (19961, any Jones matrix J 
can be transformed to a corresponding Mueller matrix M 
by use of the Kronecker product, 

M = S“^ (J (g) J*)S (4) 

where S and are the conversion matrices to transform 
between Stokes parameters and the brightness coherency 
vector. For reference they are presented below: 


S = 


0 0 \ 


/I 
0 
0 

\l -1 0 0 / 
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( 5 ) 
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( 6 ) 


Using equation 1^ the feed-error matrix D can be con¬ 
verted to a Mueller form Dm (Eq. . Each element of which 
can be directly computed using Mij = | ti{ai Jaj F), 
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/I + + |dq_>p|2) 


Dm = 


2 VI“'P—>'3 1^ T juij-s-pl 
|(|dp->qP — |dq_>pp) 


\ lra[dp^q + dp^p] 


1 “ |(Mp-><3p “ Mg-i-pP) 

R6[dp_>g -j- dq_>p] 


Im[dp^5 + di 

R,6[dp —dq —^p] Iin[dp— dq — ^p] 


1 - Re[dp^qd*^p] 
lm[dp^qd*^p] 


Im[dj 

1 + Re[dp^qd^^p 


u,p—^qdq_^p\ 


(7) 


]J 


where ai is the Pauli matrix. In this form, the {i,jY^ 
element can be understood as the response of the P*' output 
to the input. 


1.2 Polarimeter performance metric 

The polarization leakage of a dual-feed receiver is quanti- 


hed by using cross-polarization ratio (XPR) metrics (IEEE 
19981. In Carozzi & Woan (20111, a new XPR, the polarime¬ 


ter intrinsic cross-polarization ratio (IXR), was introduced. 
XPRs are used as metrics for radio communication feeds 
where the polarization of both the source and receiver is 
known. Thus an XPR can vary by choice of coordinate sys¬ 
tem ( [Carozzi fc Woan||2011[ ). The IXR is an XPR which is 
invariant under coordinate transform. This makes the IXR 
ideal for a radio astronomical polarimeter as there is no pre¬ 
ferred sky coordinate frame. The IXR in Jones formalism is 
defined as 


IXRj = ( ffmin 

ffmax 5 'min 


( 8 ) 


where ffmax and gmin are the maximum and minimum am¬ 
plitude gains of the polarimeter when using singular value 
decomposition (SVD) to decompose the system Jones ma¬ 
trix Jays. The SVD theorem (Eq. states that any Jones 
matrix J can be decomposed into two unitary transforms 
U,Vt and one diagonal transform matrix S. 


J = usv 


t _ 


u 


0 


0 


VT 


( 9 ) 


Given a noise input signal e, the sensitivity to change in 
V from Equation[^is measured by the matrix condition num¬ 
ber cond 2 (J) = k(J) = Znuix^ where Umax and (Jmin are the 

^min 

maximum and minimum singular values. An ill-conditioned 
matrix, one with a condition number much larger than 1, 
will cause an increase in the error of v with respect to the 
error of e. Conversely, a well-conditioned matrix, one with 
a condition number close to 1, will transform e into v with 
minimal effect on the error. 

Carozzi & Woan ( 2011| show that by setting the max¬ 


imum and minimum amplitude gains to be equal to the 
maximum and minimum singular values (nmax = ffmax and 
(Tmin = 5min), there is always an orthonormal choice of co¬ 
ordinates systems for the sky and the channels that gives J' 
from J such that the feed error matrix takes the form 


: + 5n 


1 1/VIXR7\ 

i/Vmu 1 ) 


( 10 ) 


The IXRj is in units of power and Equation has com¬ 
ponents of -yiXRj as a Jones matrix acts as an operation 
on an electric field vector. Intrinsic polarization leakage can 
be seen as differential gains or ‘canonical’ polarization leak¬ 
age depending on the basis. The condition number can thus 
be related back to the intrinsic polarization leakage. That 
= with a normalization factor of 


. Redehning the IXRj in terms of the condition num- 
ber. Equation [^becomes 


IXRj = 


K.(J) -(- 1 
K.(J) — 1 


( 11 ) 


This is a crucial quantity which represents a fundamental 
limit in our ability to recover the true signal. This limit is 
independent of the polarization calibration. 

The IXR is conceptually equivalent to the polconversion 
( Hamaker|2000 1 and Lorentz boost ( Britton|2000 1 transfor¬ 
mations that have been employed in previous works based 
on quaternions and geometric algebra, respectively. For ex¬ 
ample, where j3 is the velocity parameter that describes a 


Hermitian Jones matrix in Equation 11 of Britton (20001, 
k(J) = and IXRj = coth^((3). These equations enable 
meaningful comparisons between the results presented in 
this work and the notation employed in some previous stud¬ 


ies (e.g. van Straten (20131). 


To understand how the condition number, and by ex¬ 
tension the IXR, affects an observation, we can look at how 
the true sky vector e is determined. To obtain an estimate 
of the true Jones sky vector e from the observed Jones vec¬ 
tor V, the system Jones matrix Jays must be determined and 
inverted (Eq. |12|) via calibration 


I — 1 

•^sys'^ 


( 12 ) 


where e is an estimate of e due to multiple compounding 
effects. First, in the measurement of v there is a limit in 
precision due to noise. Second, Jays is not perfectly known, 
but is an estimate based on modeling and calibration. Fi¬ 
nally, the condition of the components of Jays determine how 
the errors in measurement affect the estimation of e. For an 
ill-conditioned Jays, a small error in v will result in a large 
error in estimated sky vector e compared to the true sky 
vector e. An ill-conditioned Jays matrix will lead to a noisy 
estimate of e, no matter how well known Jaya is, due to the 
inherent noise in the measurement of u. As the conditioning 
of the Jaya improves, so too does e more accurately describe 
e. 

By definition k(J) ^ 1. Ideally there is no intrin¬ 
sic polarization leakage between feeds, i.e. the matrix is 
perfectly conditioned (/max = (?min ^ k(J) = 1 and the 
IXRj —>■ 00 . That is, the two receptors are completely or¬ 
thonormal. If there is intrinsic leakage between the two 
feeds then (/max > s^min => k(J) increases and the IXRj 
decreases. In the worst case (e.g. where the two feeds are 
perfectly coupled, or one receptor’s sensitivity goes to 0), 
then (/min —>■ 0 k(J) —>■ 00 and the IXRj —>■ 1. Since the 
IXRj is a measure of feed response, it is common to use 
decibel (dB) units, IXRj^dB = 10log(IXRj). 

We can also consider the IXR in terms of Mueller ma¬ 


trices. Carozzi & Woan (20111 connects the IXR to Mueller 
matrices by showing «:(M) = «i^(J)- This relation is used to 
show the IXR in Mueller formalism is 

k(M) + 1 _ 1 + IXRj 


IXRm = 


k(M)-1 2VIXR7 


( 13 ) 
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270° 


0 ° 



IXRj (dB) 


Figure 1. IXRj of a simple ‘all-sky’ element based on the Pre¬ 
cision Array for Probing the Epoch of Reionization (PAPER) 
dipole beam model at 130 MHz. The orthogonal dipoles are ori¬ 
ented along the (45°, 225°) and (135°, 315°) axes. 


which provides a useful metric for measuring the impact of 
instrumental polarization on the Stokes parameters, espe¬ 
cially in the case of impure transformations with no corre¬ 
sponding Jones matrix. 

An example of the variation in IXRj across the field of 
view of a simple dipole element is shown in Figure The 
IXRj is maximized (70 dB) in the direction of zenith, but 
rises when observing away from boresight. The low IXRj 
structure ((0°,180°) and (90°, 270°) axes) is along the 45° 
line between the two orthogonal receptors. The variation 
in IXRj across the field of view also depends on observing 
frequency. 


Carozzi et al. (20091 and Sutinjo & Hall (20131 show the 


IXRj,dB can vary between 0 dB and 66 dB across an aper¬ 
ture array depending on pointing direction and element de¬ 
sign. For an idealized short dipole the IXRj varies smoothly 
over the observable hemisphere. But for many feeds — such 
as Vivaldi-type, bow-tie, and narrow-band half-wavelength 
dipoles — sharp intrinsic polarization leakage structures 
form across the hemisphere. Figures 3 and 5 in |Carozzi et al.| 
||2009|) show the IXRj over a hemisphere for short dipoles 
and Vivaldi-type elements. Figures 2 and 3 in [Sutinjo "fcl 


Hall (20131 show the IXRj across the held of view (FoV) 


of an MWA bow-tie element. These published values and 
maps provide insight into what range of IXRj to use in our 
simulations. 


1.3 IXR and signal-to-noise ratio 

The error on a ToA measurement is, in general, a function of 
the S/N of a given observation. By S/N, we hereby refer to 



IXR (dB) 

Figure 2. The average observed S/N of J1603—7202 as a fraction 
of the ideal S/N for a range of IXR values using ‘gain’ (blue/solid) 
and ‘full’ (red/dashed) calibration methods described in Section 

[2:21 


the peak pulse value in the Stokes I prohle to the standard 
deviation of the off-pulse signal. Under ideal circumstances, 
the S/N increases with the square root of integration time. 
For a polarized source, instrumental intrinsic polarization 
leakage will result in a lower observed S/N compared to the 
ideal S/N, i.e. that obtained from a system with no intrin¬ 
sic polarization leakage, for a given amount of integration. 
The blue/solid line in Figure|^shows the fractional observed 
S/N compared to the ideal S/N as a function of the IXR for 
J1603—7202. The red/dashed line shows how using the in¬ 
verse of a poorly conditioned matrix for calibration amplifies 
the noise in the measured profile. This will be discussed fur¬ 
ther in Section All the simulated pulsars have a similar 
response. An effect of a low IXR is the introduction of a 
differential gain between feeds. As the IXR goes to 0 the re¬ 
ceiver system becomes effectively blind to one polarization, 
thus the observed S/N is approximately half (blue/solid line) 
that of the ideal S/N in the limit IXR —>■ 0. When calibra¬ 
tion is applied, not only is there a differential gain effect, but 
the inversion of the ill-conditioned matrix will significantly 
degrade the S/N of any profile (red/dashed line). The ideal 
S/N is achieved only in the limit as IXR —>■ 00 , and there 
is a one-to-one correspondence with integration time. A ref¬ 
erence integration time of Tint = 1 is defined as the time it 
would take to build up an ideal S/N of 1000. All integration 
time values quoted in this paper are a fraction of this ref¬ 
erence integration time. The relationship between Tint and 
ideal S/N SNi is 



2 THE SIMULATIONS 

We have performed simulations with the goal to quantify 
the effect of intrinsic polarization leakage on pulsar ToA 
measurements. We have sampled the three-dimensional pa¬ 
rameter space that includes the IXRj, ideal pulse profile 
S/N, and calibration error, which covers current telescope 
measurements and future telescope specifications. For each 
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point in the sampled parameter space, ‘observed’ profiles 
are generated for 500 epochs, by modifying a template pro¬ 
file with the appropriate Mueller matrices. For every epoch, 
we stochastically generate the Jsys Jones matrix, and the 
corresponding Mueller matrix Mays, for a given intrinsic po¬ 
larization leakage. The form of this matrix is described in 
§2.2[ The observed profile is then calibrated by multiplying 
by the inverse of the system Mueller matrix with additive 
random calibration errors. 

A ToA is determined at each of the 500 epochs, using a 
standard timing method (section 2.31 included in PSRCHIVE 
(Hotan et al. 2004 van Straten et al. 20121. In a normal 
pulsar timing experiment, a model would then be fit to the 
ToA measurements using TEMP02 (Hobbs et al. 20061 and, 
the goodness of the model would be measured by the root 
mean square (rms) of the timing residuals. Since we are us¬ 
ing a simple model of an isolated, stable pulsar of constant 
period, this rms represents the ideal rms for a given set of 
parameters. 

In the following we give further details of the steps of 
the simulation. 


2.1 Pulse profiles 


We performed our simulations using profiles based on the 
mean pulse profiles of 20 well studied MSPs from the Parkes 


Pulsar Timing Array Manchester et al. (20131. It is not our 


intention to reproduce the results of that work. We use these 
MSPs as they cover a wide range of profile and polariza¬ 
tion structures, and perform our simulations to much higher 
S/Ns than is possible with current experiments. 

A high S/N version of each reference profile was gener¬ 
ated by applying a rolling Hann filter with a width of ~ 1% 
of the profile. This has the effect of removing high-frequency 
components and reduces the amplitude (up to a few percent) 
of narrow profile structures. This was done in order to avoid 
very low uncertainties in our TOA measurements due to the 
presence of the same high frequency structures in the tem¬ 
plates and the simulated noisy data. Though we retain the 
pulsar names, by applying this filter, the new profiles are 
approximations to the original profiles and are meant only 
to be a representative set of profile structures. The low fre¬ 
quency components that remain, may cause our uncertain¬ 
ties to be systematically lower than real timing experiments 
for these pulsars, however this effect is consistent across our 
simulations and therefore does not affect our final result. 

The ideal profiles are then used to generate imperfect 
profiles at each epoch using the process described in Sec¬ 
tion [2?^ These profiles were also used as the template pro¬ 
files when performing timing. The first and third columns 
of Figures at the end of this paper are the Stokes 

parameters of the ideal profiles. 


2.2 Simulated observations 

Our simulated observations produce an estimated sky Stokes 
vector e® for a given pulsar per epoch by estimating the 
system Mueller matrix, including calibration error, intrinsic 
polarization leakage, and S/N parameters. Figure shows 
the stages to arrive at an example e®. 

The explicit RIMEs defined in Equation can be sim¬ 
plified for the simulations. By using the mean profile we 


are making the simplification that the bandpass B and time 
varying G system gains have been solved for and applied to 
the observed signal. Thus, B and G are unity diagonal matri¬ 
ces. We include a polarization calibration error term AJ into 
our system to simulate the effect of imperfect calibration of 
B and G. The nominal feed matrix C is a telescope-specific 
basis transform, and will not affect the IXR as it is a co¬ 
ordinate independent metric (i.e. IXRc = oo). Thus, the 
polarization leakage matrix D is the only matrix which will 
vary in our parameter space. In practice D is a function of 
frequency, but as the profile is an average across a frequency 
band, so too is the IXR a frequency averaged intrinsic po¬ 
larization leakage. 

To generate a simulated observed profile we start by cre¬ 
ating a intrinsic polarization leakage Jones matrix represen¬ 
tation. We define the polarization leakage to be 1/VlXRj, in 
decibel units the intrinsic polarization leakage is —IXRj,dB. 
We have chosen this definition as polarization leakage is a 
common concept within the community. The IXR is a mea¬ 
sure of both the cross-coupling between receiver feeds, which 
is typically thought of as ‘polarization leakage’, and the ap¬ 
parent leakage due to differential feed gains and thus is a 
more complete metric for ‘polarization leakage’ over previ¬ 
ous definitions. A higher intrinsic polarization leakage im¬ 
plies the two feeds are more coupled together than a lower 
intrinsic polarization leakage. To sample a broad range of in¬ 
trinsic polarization leakage values we sample the space 0 dB 
to — 30dB, where the upper limit of 0 dB is effectively blind 
to one polarization, such as a single polarization receiver. 
The intrinsic polarization leakage sample space in IXR no- 


the condition number is related to the IXRj is 


tation is 1 ^ 1/vTXRj ^ 0.031. By inverting equation 


11 


k(J) 


yrxRj + i 

VIXR7- 1 


(14) 


For each run of the simulation with a given set of param¬ 
eters, we construct a system Jones matrix J^yg (Eq. i by 
generating a random complex matrix from a normal distri¬ 
bution (^ = 0,(T = 1) for D. All other matrices in Jays are 
simplified to unity matrices as discussed earlier in the sec¬ 
tion. The random matrix is decomposed using SVD (Eq. 
[^, and the singular values in S are replaced with those for 
the simulated IXR parameter. Without loss of generality we 
normalize S using CTmax = (?max = 1. The normalized condi¬ 
tion number becomes k(J) = l/gmin and the system Jones 
matrix due to intrinsic polarization leakage becomes 


Jays = USVt = U vt (15) 

We use a random matrix as there is an infinite set of Jones 
matrices for a given IXR. The decomposition by SVD, and 
reconstruction steps are to maintain the same scaling as with 
the calibration error Jones matrix. 

For error in the polarization calibration, we use a sample 
space starting at the ideal perfect calibration (0%) up to 
15% calibration error. This calibration error represents the 
imperfect measurement of B and G, the inverse of which 
are applied to the observed signal to integrate the profile 
in both time and frequency. A calibration error matrix AJ 
is generated from a random normal complex distribution 
{n = 0,<y — precent error/100). The estimated system Jones 
matrix is Jays = Jays +AJ. This error parameter space covers 
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Figure 3. Stages of simulation for J1603—7202. Central third 
of pulse profile plotted. Stokes Parameters: I (black/solid), Q 
(green/dashed), U (red/dotted), and V (cya n/solid) |(a)] True sky 
Stokes vector of the pulsar profile, (b) True system Mueller 
matrix Maya with the condition number k(J) set for IXRj ub = 7 
dB applied to e’®. (c) Added receiver noise for an SNi of 100. (d) 
Estimated Stokes vector e"® of the pulse profile used in timing 
exp eriment using ‘gain’ calibration with 20% calibration error. 


(e) Estimated Stokes vector of the pulse profile used in timing 


experiment using ‘full’ calibration with 20% calibration error. 


a wide range of polarization calibration situations, although 
the typical calibration error is of the order of a few percent 
(Han et al. 20091. In practice, the polarization calibration 
solution will also be affected by the IXR. We simulate this 
effect by using two different calibration techniques discussed 
later in this section. Additionally, this error does not include 
the potential error from approximating D based on beam 
modelling and observations that would be done in a real 
system. Using equation Jsys and Jays are converted to 
their Mueller forms Mays and Maya- 

To set the simulation S/N we add a system noise com¬ 
ponent appropriate for the integration time of the simulated 
observation as described in o This is set by scaling the 
ideal profile e®, which has been normalized such that the 


Stokes / peak is unity, by a scalar SNi value and adding a 
system noise Stokes vector n® which is a set of real random 
values from a normal distribution (/i = 0, cr = 1). The noise 
component n® is a direction-independent effect. This is a 
practical approximation when the system noise is receiver 
noise dominated or when the sky noise is isotropic on the 
scale of the beam primary lobe. In the extreme case where 
all the system noise is direction-dependent in the direction 
of the source the effect on rms timing will only be a loss in 
S/N as seen in Figure]^ As we will see in the following sec¬ 
tions, the effect of direction-independent noise is to further 
degrade the timing solutions beyond the loss in S/N. For 
simplicity we are using only direction-independent noise for 
the simulation. 

An ideal S/N is set to be in the range vTiF to VTCF 
for our simulations, though we note: the observed S/N is 
reduced as the condition number increases, and thus the 
actual S/N is a function of IXR as shown in Figure]^ The 
estimated Stokes vector e® of the observed pulsar is 


e® = M. 


(sNi 


X Ma 


e® + n® 


(16) 


The inversion of Mays is the calibration stage (Eq. 121. 


We have simulated two types of calibration. The first is 
a ‘gain’ calibration where we are only interested in solv¬ 
ing the bandpass and electronic gain of the system, that is 
Jays = BGC + AJ = I -I- AJ, where I is the identity matrix. 
This is an ideal calibration solution where the gain terms can 
be solved for independently of any D effects, such as if using 
a known noise reference in a single dish system. The second 
type of calibration is a ‘full’ calibration where the D term 
is included. Jays = BGCD-b AJ = D -|- AJ. This is a more 
realistic approach, as the IXR will affect any calibration so¬ 
lution. And, when performing timing with an array, complex 
gain solutions are required to combine signals in phase by 
using a sky calibrator source or self-calibration. The cali- 
bratability of an array is a topic which should be studied in 
further work. Figurej^shows the effect of these two methods 
on observed S/N, independent of pulsar and ideal S/N. Ex¬ 
amples of these type s of c alibr ation on the observed profile 
are shown in Figures 3(d) and 3(e) After a simulated e® is 


produced, a ToA is then determined with standard pulsar 
timing software. The effect of these calibration methods on 
timing will be shown in the followingsection. The simulation 
code is available as a git repositorjIM 


2.3 Methods for TOA determination 


For ToA measurements, three standard methods are used: 
total intensity ( |Taylor||1992[ ), invaria nt interval ||Britton| 
20001, and matrix template matching (van Straten 20061. 


These methods are included in PAT, from the PSRCHIVE pack¬ 
age. For the total intensity method, each observed Stokes / 
profile is cross-correlated with the ideal template. The in¬ 
variant interval technique uses all Stokes parameters in the 
form of a Lorentz 4-vector {F — . By in¬ 

cluding all Stokes parameters complete information is used. 
However using the invariant form, the S/N decreases, lead¬ 
ing to a less precise ToA determination compared to the 


^ https://github.com/grifHnfoster/pulsar-polarization-sims 
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Figure 4. Contour plots of time of arrival rms noise (fis) from simulation of PSR J1603—7202 as a function of calibration error and 
intrinsic polarization leakage. The top row shows the effect of intrinsic polarization leakage on timing for the three methods when ‘gain’ 
calibration is applied. The bottom row is the same simulations but with ‘full’ calibration applied. There is a strong dependence on 
polarization calibration error when using the total intensity method. The solutions when using the invariant interval are effectively 
identical for both types of calibration. Matrix template matching outperforms the invariant interval method in both the calibrated and 
uncalibrated cases. An ideal SNR of 100 was used for these figures. 


intensity fitting, when the source is highly polarized or in¬ 
trinsically weak (see [van Straten| ( |2013[ ) for more informa¬ 
tion on this effect). Matrix template matching represents 
the profile in Jones notation and produces a ToA measure¬ 
ment while simultaneously solving for a calibration solution 
by transforming the template profile to the observed profile. 
In practice, the S/N requirements for timing often “force” 
a pre-processing step of frequency and time averaging. This 
requires some prior level of polarization calibration, as de¬ 
scribed in the previous section, rendering the assumptions 
for matrix template matching no longer strictly valid, by 
introducing covariances between the Stokes parameters, es¬ 
pecially fir large polarization leakage. 

Figurej^shows the ToA rms derived from the three tim¬ 
ing methods used in simulations of the well studied MSP 
J1603—7202. There is a strong dependence on polariza¬ 
tion calibration error when using the total intensity timing 
method. Applying the ‘full’ calibration results in better tim¬ 
ing solutions except in the high intrinsic polarization leak¬ 
age region (> —5 dB), compared to using only the ‘gain’ 
calibration. We will use the ‘full’ calibration simulations for 
our total intensity results. The timing error, when using the 
invariant interval method, is effectively the same for both 
types of calibration. 

For the matrix template matching method there is a 
significant change in the ToA rms at high intrinsic polar¬ 
ization leakage between the two types of calibration. This 
is due to ‘full’ calibration introducing covariances between 
Stokes parameters, compared to just ‘gain’ calibration. Sim¬ 
ulations using ‘gain’ calibration and timed with matrix tem¬ 
plate matching produce the best timing results for a given 
calibration error and intrinsic polarization leakage. We will 


present matrix template matching results using both types 
of calibration, but will focus our analysis on the ‘full’ cal¬ 
ibration as is represents a more pragmatic result in terms 
of system calibration, especially in when using an array for 
timing which requires multiple signals to be combined in 
phase. 

In practice, total intensity and invariant interval are 
commonly used methods, where as the matrix template 
matching method is more rarely used due to the difficul¬ 
ties in practical polarization calibration. However, as timing 
limits are pushed, matrix template matching will become 
necessary Ivan Straten|2013 1. 


3 RESULTS 

As polarization calibration error is not the focus of these 
simulations, we can compress the contour plots in Figure]^ 
into a more concise and useful figure by collapsing the cali¬ 
bration error axis. Figure shows the bound range for each 
method on a single plot. The narrowness of the bounds for 
the invariant interval and matrix template matching show 
their independence from polarization calibration error. Plots 
of this style for all the simulated MSPs are shown at the end 
in Figures - 110| These figures show that as the intrinsic 
polarization leakage decrease, so too do the timing residuals 
using all methods. 

We have chosen to use the ToA rms noise as the metric 
for these results. This is a measure of the time of arrival scat¬ 
ter for the simulations. Using the average ToA uncertainty 
ctoA underestimates the error. As per Eq. 8.2 of |Lorimer| 
| |2005[ ) cttoA — W/SNo, where W is the pulse width and 
SNo the observed S/N. In our simulation results we see the 
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Figure 5. Time of arrival rms noise for J1603—7202 as a 
function of intrinsic polarization leakage using the three tim¬ 
ing methods: total intensity (red/no border), invariant interval 
(gray/dashed border), matrix template matching (full calibration: 
orange/dotted border, gain calibration: blue/dot dash border). 
The width of the lines show the effect polarization calibration 
error has on the ToA rms noise, a polarization calibration error 
from 0% to 15% was used. These simulations use a fractional inte¬ 
gration time of 0.01 which would produce an ideal signal-to-noise 
ratio of 100. 


Figure 6. Time of arrival rms noise using the total inten¬ 
sity method as a function of intrinsic polarization leakage for 
J1603—7202. Lines cover a range of polarization calibration er¬ 
rors from 0% (blue) to 27% (red). A normalized integration time 
of 0.01 was used for all the lines. A counter-intuitive effect occurs 
when the intrinsic polarization leakage is high, where the rms 
noise becomes inversely proportional to the calibration error be¬ 
cause calibration errors improve the condition of the calibration 
matrix Msys. 


intrinsic polarization leakage affects timing results beyond a 
simple reduction of the S/N, producing poorer rms timing 
residuals compared to the average cttoA of the same sim¬ 
ulated observations. In the limit the IXR,j —>■ oo, the rms 
noise will reach utoA- 

We computed the reduced solution to measure the 
goodness of fit for each sample of the parameter space. When 
computing the of th® measured time of arrival against 
the expected time of arrival, the low intrinsic polarization 
leakage case results in a good fit, as would be expected. As 
the intrinsic polarization leakage increases the fit degrades. 
Paradoxically, as the intrinsic polarization leakage increases 
to the highest values, the reduced fit approaches 1. This 
is because cttoA grows exponentially large, and the error of 
the fit is within the limits of the ToA variance. 

Although the diverse set of profiles we have simulated 
produce different results, seen in Figures - 110| there are 
general trends we observe. All the simulations with the ‘full’ 
calibration profiles show exponentially increasing rms tim¬ 
ing solutions as the intrinsic polarization leakage increases. 
Using only ‘gain’ calibration produces good timing results 
in this region, but is not as realistic as the ‘full’ calibration 
simulations. 

Generally, using matrix template matching, with ‘gain’ 
or ‘full’ calibration, produces better timing solutions com¬ 
pared to other methods at all intrinsic polarization leakage 
values. For the majority of the ‘gain’ calibrated profiles, ma¬ 
trix template matching produces timing solutions that are 
only weakly dependent on intrinsic polarization leakage. It 
is worth noting that, for a few profiles, such as J1022-I-1001, 
matrix template matching of the ‘full’ calibration profile pro¬ 
duces lower rms timing residuals than the ‘gain’ calibrated 
profile. This profile has highly polarized components, which 
are distorted by high intrinsic polarization leakage. At high 
intrinsic polarization leakage, the profiles with ‘full’ calibra¬ 


tion result in high rms timing residuals and dependence 
on polarization calibration error when using matrix tem¬ 
plate matching. However, at low intrinsic polarization leak¬ 
age both types of calibration converge to a similar timing 
solution. 

The results from using the invariant interval method 
are consistent across all profiles. In a few cases, for example 
J1744—1134 and J1824—2452, the invariant interval method 
significantly underperforms compared to the other methods. 
Both these pulsars are highly polarized; therefore, in the 
invariant interval, the power in Stokes I is almost entirely 
cancelled out by the other Stokes parameters. 

Timing solutions with the total intensity method have 
a notable dependence on the polarization calibration error. 
Plotting the bound regions as individual lines for different 
polarization calibration errors in Figure there is a polar¬ 
ization leakage point (e.g. for PSR J1603—7202 it is around 
—4 dB), after which systems with higher polarization cali¬ 
bration error produce better timing solutions than an ideally 
calibrated system. This is true for all MSP profiles used in 
our simulation, though transition points vary between —3 
dB and —7 dB. In this region of the parameter space, this 
counter-intuitive effect comes about because the error-free 
calibrator transformation also has high intrinsic leakage and 
therefore further reduces the observed S/N of the calibrated 
profile. Adding polarization calibration error reduces the in¬ 
trinsic polarization leakage, leading to an increase in the 
observed S/N. There is potential use for the total intensity 
method at high intrinsic polarization leakage when only us¬ 
ing a ‘gain’ calibration. Needless to say, this regime should 
be avoided and this sets a maximum limit on the allowable 
intrinsic polarization leakage to around —5 dB. 

Figure shows the residual timing rms of J1603—7202 
as a function of S/N and IXRj using a calibration error 
of 5%. For a given integration time, the achievable time of 
arrival rms noise is dependent on the system polarization 
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Figure 7. Contour plot of time of arrival rms noise (/rs) using 
matrix template matching as a function of intrinsic polarization 
leakage and integration time for PSR J1603—7202. Integration 
time is a fraction of the reference integration time, see 



Figure 8. Time of arrival rms noise of PSR J1603—7202, using 
the matrix template matching method, as a function of integra¬ 
tion lengths for a range of intrinsic polarization leakage values. 
Horizontal lines indicate rms thresholds: 0.5 ps (dotted), 1 ps (dot 
dash), 5 ^s (dashed). 


value, e.g. in simulation of J1603—7202 for Tint = 0.1 a sys¬ 
tem with —25 dB intrinsic polarization leakage can achieve a 
ToA rms in the timing residuals of around 300 ns compared 
to 3000 ns with a —5 dB intrinsic polarization leakage sys¬ 
tem. This variation due to intrinsic polarization leakage can 
be better seen in Figure]^ which shows the residual rms as a 
function of integration length, for a range of intrinsic polar¬ 
ization leakage values. The rms will continue to decrease as 
the integration time increases. The system intrinsic polar¬ 
ization leakage sets the required integration time to achieve 
a desired rms. In we discuss the importance of Figure 
in optimizing science return when there is limited available 
observing time. 


4 DISCUSSION 

Gravitational wave detection using high precision timing of 
MSPs constitutes a key science project for the SKA. This 
imposes a requirement on the polarization specifications. In 
|Cordes et al.| | |2004[ ), a case is presented which sets the re¬ 
quired polarization purity level to —40 dB to accomplish the 


SKA pulsar key science goals. This polarization purity value 
is different from what we have defined as polarization purity; 
it is a measure of the final, calibrated Stokes data and not a 
specification of the front-end design as we have considered 
here. 

Figure shows that for a given intrinsic polarization 
leakage, a desired timing residual rms can be achieved with 
sufficient observing time. This, of course, ignores the other 
systematic effects that are part of a timing observation and 
does not include the additional sensitivity modulation of 
the primary beam shape. We have only focused on intrinsic 
polarization leakage, which is an effect on any feed design. 
The main issue is that the intrinsic polarization leakage has a 
strong effect on the required observation time, which is a lim¬ 
ited commodity. For our simulation of MSP J1603—7202 in 
Figure|^the difference in required observing time to achieve 
a desired rms noise at —15 dB compared to —30 dB is a fac¬ 
tor of 1.5. With limited available observing time, we would 
like to set an upper limit on the intrinsic polarization leak¬ 
age above which it is no longer optimal to be using obser¬ 
vation time on a measurement. Carozzi et al. (20091 and 
|Sutinjo fc Hall| ( |2013[ ) show the IXRj^b for typical feeds to 
be somewhere between 30 dB and 66 dB at boresight. We 
see that observing with an intrinsic polarization leakage of 
— 15 dB implies at least a 50% increase in observing time 
compared to that of a typical feed. A high intrinsic polar¬ 
ization leakage is not a design issue for a ‘classic’ single pixel 
dish system, where a low intrinsic polarization leakage can 
be achieved when observing a source on axis. This is not the 
case with aperture arrays, PAFs, and multi-beaming with 
single pixel dishes. In these cases, the source will likely not 
be located in the optimal intrinsic polarization leakage re¬ 
gion of the beam. For aperture arrays, a source will rarely, if 
ever, be on axis. Returning to the example beams in|Caroz^ 


et al. (20091 and Sutinjo & Hall (20131, we see that the po¬ 


larization leakage values can quickly increase to above —10 
dB away from zenith. This effectively limits the declination 
range of sources, depending on the array latitude, for pulsar 
timing. 

From our simulation we see there is an intrinsic polariza¬ 
tion leakage lower limit on feed design at which point there is 
minimal return in terms of reducing the timing residual rms 
with further engineering investment, for a given integration 
time. As there is a cost to every incremental improvement in 
IXR, we would like to present our simulation results in terms 
scientific return for marginal improvements in engineering 
specifications. In an effort to create a meaningful engineer¬ 
ing intrinsic polarization leakage lower limit. Table lists 
the fractional improvement in ToA rms noise for different 
IXR values. We have picked a typical S/N (rint = 0.01) for 
a timing observation. Columns 2, 4, 6, and 8 are the timing 
residual rms for each MSP at IXRj = 10, 20,30,40 dB re¬ 
spectively. Columns 3, 5, and 7 are the percentage change 
in the rms with the changes in IXR. This table shows the 
diminishing marginal utility of improving IXR for the ben¬ 
efit of decreasing the time of arrival rms noise. There is, on 
average, a 29% improvement in the rms when improving the 
IXRj from 10 dB to 20 dB, but can vary significantly with 
profile shape. For example, timing of J0711-6830 is largely 
uneffected by improvements in IXR, while timing of J1603- 
7202 improves with each increase in IXR. Going from 20 dB 
to 30 dB there is a small improvement, but going above 30 
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dB provides essentially no improvement. This indicates that 
for a feed with IXRj,dB > 30 dB there is limited fractional 
improvement in pulsar timing capabilities. It may be worth 
considering that low intrinsic polarization leakage across the 
held of view may be preferable to optimizing for minimal in¬ 
trinsic polarization leakage around boresight. 


5 CONCLUSION 

On the pathway towards the SKA a number of aperture 
arrays, dishes, and PAFs are being developed as precursor 
instruments. As pulsar timing is a key science project, design 
of these instruments should take into account the intrinsic 
polarization leakage specihcation we have presented in this 
paper. 

There is a relative increase in required integration time 
as a function of the feed IXR, as seen in Figure]^ At high 
intrinsic polarization leakage this can make the required in¬ 
tegration time inefHciently long. 

We have shown that there are diminishing returns (Ta¬ 
ble on building feed systems which have intrinsic polar¬ 
ization leakage below —30 dB in the direction of observation. 
Achieving this intrinsic polarization leakage limit should be 
easily affordable for single pixel dishes on axis, where the 
leakage is at a minimum. However, aperture arrays, PAFs, 
and multi-beam systems, where observations are not always 
made in the direction of minimum leakage, could have lim¬ 
ited use for pulsar timing if intrinsic polarization leakage 
is not taken into account while developing the feed system. 
The complex aperture array and PAF beams with sharp, fre¬ 
quency dependent structure lead to regions of high intrinsic 
polarization leakage. 

Given the effect of intrinsic polarization leakage on pul¬ 
sar timing and a costing model for a feed design, a desirable 
optimization could be to maximize the average IXR across 
the intended usable field of view for the element and not 
just in the direction of boresight. 

The calibratability of an array has a key effect on tim¬ 
ing as we have shown with the idealized ‘gain’ calibration 
technique against the more realistic ‘full’ calibration. Be¬ 
yond this work, there is scope for additional work on effect 
of calibration on timing. Additionally, the matrix template 
matching method should be extended to account for the co- 
variances between the Stokes parameters induced by poorly 
conditioned calibration matrices. 

Ideally, we can further refine these values on precursor 
instruments. The MeerKAT, KAT-7, ASKAP, and APER- 
TIF arrays will provide a platform to study dish array polar¬ 
ization effects with multi-beam feeds and PAFs. LOFAR and 
MWA, though not ideal for pulsar timing experiments due 
to the low observing frequencies, will be useful to study in¬ 
strumental polarization in aperture arrays. A study of these 
array polarization properties is necessary to assure the SKA 
science specifications can be met for pulsar timing. 
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Figure 9. Time of arrival rms noise from simulation of various MSPs. The first and third columns are the Stokes parameters of the 
ideal profile: I (black/solid), Q (green/dashed), U (red/dotted), and V (cyan/solid). The second and fourth columns show ToA rms noise 
(ps) as a function of intrinsic polarization leakage when using total intensity (red/no border), invariant interval (gray/dashed border), 
and matrix template matching (full calibration: orange/dotted border, gain calibration: blue/dot-dash border) methods on the profile. 
The width of the lines show the effect polarization calibration error has on the rms, a polarization calibration error from 0% to 15% was 
used. These simulations use a fractional integration time of 0.01 which would produce an ideal signal-to-noise ratio of 100. Plots continue 
in Figure [To| 
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Table 1. Fractional improvement in time of arrival rms noise between two IXRj values for the simulated MSPs using the matrix 
template method. Percent change is computed as AtoA = 100(rmsi — rmsf)/rmsi, where rmsi is the initial (lower) IXRj ToA rms noise, 
and rmsf is the final (higher) IXRj ToA rms noise. These simulations used a fractional integration length of 0.01, which is an ideal SNR 
of 100. The reported rms values have an uncertainty of ~ 1% set by the number of ToA simulations (n = 5000). The last rows of the 
table are the average ToA rms noise percent improvement and minimum/maximum range. 
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Figure 10. Continuation of Figuresee that figure for description. 






































































































